Cite this Article
Shao Ming-Zhe, Wang Yan-Ting, Zhou Xin. Fast and accurate determination of phase transition temperature via individual generalized canonical ensemble simulation. Chinese Physics B, 2020, 29(8): 080505  
Permissions
Fast and accurate determination of phase transition temperature via individual generalized canonical ensemble simulation
Shao Ming-Zhe1, Wang Yan-Ting2, Zhou Xin3, †
College of Light Industry Science and Engineering, Tianjin University of Science and Technology, Tianjin 300457, China
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China

 

† Corresponding author. E-mail: xzhou@ucas.ac.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11574310, 11674345, and 21733010) and Beijing National Laboratory for Molecular Sciences, China (Grant No. BNLMS201835).

Abstract

It is very important to determine the phase transition temperature, such as the water/ice coexistence temperature in various water models, via molecular simulations. We show that a single individual direct simulation is sufficient to get the temperature with high accuracy and small computational cost based on the generalized canonical ensemble (GCE). Lennard–Jones fluids, the atomic water models, such as TIP4P/2005, TIP4P/ICE, and the mW water models are applied to illustrate the method. We start from the coexistent system of the two phases with a plane interface, then equilibrate the system under the GCE, which can stabilize the coexistence of the phases, to directly derive the phase transition temperature without sensitive dependence on the applied parameters of the GCE and the size of the simulation systems. The obtained result is in excellent agreement with that in literatures. These features make the GCE approach in determining the phase transition temperature of systems be robust, easy to use, and particularly good at working on computationally expensive systems.

PACS: ;05.70.-a;;05.70.Fh;;02.70.Ns;
Keyword:;phase transition;;enhanced sampling;;metastable state;;molecular dynamics;
1. Introduction

The first-order phase transitions, such as the water/ice transition, are very popular phenomena in nature. Molecular modeling and simulations offer microscopic insights to both the thermodynamics and kinetics of phase transitions.[1–11] The phase coexistence temperature is one of the center values to verify models, and efficiently calculating the temperature via simulations to compare with experiments is key in the molecular modeling and simulations of phase transitions.

Estimating the phase transition temperature via simulations has been a long term issue. Take the TIP4P water model for example, it was proposed by Jorgensen in 1983[12] and well performs in revealing the density of liquid water, the density near the critical point, and the enthalpy of vaporization.[13] A rough estimate of the melting point of ice Ih was first given by Kroes[14] as 230 K < Tc < 250 K for the water model by monitoring the translational mobility of atoms on the ice/water interfaces. Later in 2000, with free energy (FE) calculation, Gao indicated Tc = 238 ±7 K.[6] While another study proposed by Vega in 2005 implied a similar result, Tc = 232 ±5 K.[15] Although the FE approach determined phase transition temperature is theoretically flawless, refers to the state at which two phases have Gibbs free energy identical, tiny uncertainty in free energy estimation can lead to a totally different result. Gao described this as “1 percent in relative error leads to a shift of melting point more than 10 K”, the FE approach seems to reach its limit.

Direct coexistence (DC) simulation of phase transition is also a generally accepted method to determine the phase transition temperature. This method is particularly popular in extracting the gas–liquid coexisting line on the phase diagram,[16–19] and also can be applied in liquid/solid systems.[5,20,21] Those works of determining the phase transition temperature with DC approach show great agreement with the FE approach, and sometimes better precision if the data sample is sufficient. For example, with the water model TIP4P-2005 (by reparameterising from the TIP4P),[22] Conde calculated the melting point of the ice Ih by both FE and DC approaches, that is, 252±5 K and 249±3 K, respectively.[21] In 2017, by performing simulations covering a range of temperature space, repeating the simulation with different seeds for 5 times, using the potential energy evolution in tens of nanoseconds as the phase transition indicating order parameter, Conde developed the DC approach with exhaustive study, obtaining Tc = 249.5 ± 0.1 with the best precision to date.[5] However, the plenty of simulations and expensive CPU cost have limited the application of this approach.

The Hamiltonian Gibbs–Duhem integration offers another option. A number of studies have demonstrated that once the melting point of a model is known, the melting point of a different model can be easily estimated with the Gibbs–Duhem integration.[10,15] However, this reference model dependent estimate is not universally applicable, one can neither extract the melting point of ice Ic from an ice Ih reference, nor a silicon phase transition reference.

In this work, we aim to provide a simple and universal method with high precision and efficiency for determining the phase transition temperature. To achieve this, the generalized canonical ensemble (GCE)[23–28] has been implemented successfully in Lennard–Jones (LJ) particles and ice/water systems with the mW, TIP4P-2005, TIP4P-ICE water models. Unlike in the normal canonical ensemble simulations where phase-coexistence states are unstable thus can not be focused on, the GCE allows sufficiently sampling in the phase-coexistence regions, making a single simulation of DC approach be enough to calculate the phase transition temperature with considerable precision.

This work is organized as follows. Section 2 presents the theory for determining the phase transition temperature. Section 3 describes the models and method used in this work. Section 4 presents the results in various systems. The papers ends with a discussion and conclusions.

2. Theory
2.1. Generalized canonical ensemble

The GCE method has been explicitly demonstrated in previous studies,[23–25] here we present a brief overview of this method. In GCE, the energy-dependent thermostat temperature is applied to couple with the system, in physics, it is equivalent to apply a finite-size thermostat so that the temperature of the thermostat varies as the energy of the system (and the energy of the thermostat, since the total energy of the system and thermostat is constant). In GCE, we choose the simplest form of the energy-dependent (reverse) temperature

Here with the Boltzmann constant kB (setting as unit in this work) and a constant temperature parameter T0. The other two constant parameters in the GCE are α and E0, and we set α ≥ 0. The corresponding conformation distribution in the GCE is

equivalent to the normal canonical ensemble under an effective potential energy surface

Here U(r) is the potential energy of the original system, r is a simple denotation of the high-dimensional conformation, such as all atomic position vectors, and the integral about r is over the whole conformational space. Obviously, the GCE is an extension of the normal canonical ensemble where α = 0. Similar to the GCE with constant volume, we have the corresponding GCE with constant pressure, denoted as gNPT, where the inverse temperature is written as

with H = E + PV being the enthalpy of the system.

The GCE (or gNPT) is equivalent to applying an effective potential energy, thus its implementation is simple and direct by rescaling the physical force on each atom. It is also possible to reset the temperature of the thermostat at each MD simulation step under the original potential energy surface, from the viewpoint of applying a finite-size thermostat. The two implementations are thought to be equivalent to each other for achieving the equilibrium properties.[25,26]

The GCE aims to use a large positive α to visit sufficiently the phase-coexistent conformational regions, while the canonical ensemble (α = 0) can not focus in the regions due to the thermodynamical instability of the coexistent states. We define the microcanonical entropy S(E) = ln ∫ drδ(E – U(r)), and will show that S″(E) = ∂2 S/∂ E2 > 0 is the condition of the thermodynamical instability. The probability distribution of GCE in the energy space is

Here Em is the maximal point of the probability distribution, thus the simulation focuses on visiting the neighboring energy region of Em when and only when α > S″(Em). The Em satisfies the equation

i.e., the cross point of the curve S′(E) = ∂S/∂E which is the inverse interior (system) temperature and the curve β(E) which is the inverse exterior (thermostat) temperature. Obviously, if α < S″(Em), Em is the minimal point of the probability distribution, thus the neighbor of Em is unstable in the GCE, which can not be sufficiently visited. In the phase-coexistent regions, S″(E) > 0, thus they can not be visited sufficiently in the canonical ensembles (α = 0) but can be stabilized in the GCE by using larger α > S″(E).

If neglecting the possible small asymmetry of the GCE probability distribution around Em, we approximately have , the average energy in the GCE. Thus we can interpolate the whole curve S′(E) via the obtained points in a series of GCE simulations with different E0 (or β0). The S′(E) curve is sufficient to get thermodynamical properties of the original system in the canonical ensemble with any temperature. Actually, in the pure phases, the average energy in a GCE, , is approximately equal to that in the canonical ensemble with the corresponding temperature .

2.2. Phase transition temperature

At the phase transition temperature (or named as the phase coexistence temperature) Tc, the two pure phases have the same Helmholtz (or Gibbs) free energy as a function of temperature, As(Tc) = Al(Tc). Here we use the Landau free energy F(E; T) = E – T S(E), and . Here we neglect the varying of distribution width of energy on temperature, and is the mean energy in the canonical ensemble with the temperature T. Thus, we can determine Tc from the equation El – Tc S(El) = Es – Tc S(Es). It is nothing else but the well-known Maxwell equal area rule

Here Es and El are the cross points of the curve S′(E) with βc = 1/Tc, i.e., the mean energies of the l and s phases at Tc. It means that

Therefore, when , the system is in the s phase with the mean energy Es, but when , it is in the l phase with the mean energy El. When , the two phases have the same free energy, and the states with the energy in [Es, El] are metastable (superheating and supercooling) or unstable (phase coexistence) in the canonical ensembles. It is worth to mention here that E should be replaced by the enthalpy H if using the ensemble with the constant pressure instead of the constant volume.

Usually, ones can run simulations in canonical ensembles to estimate the temperature Tm at which the system often melts from solid to liquid phase, and the temperature Tf at which the system easily freezes from liquid to solid. The coexistence temperature Tc is between the two temperatures, and is possibly estimated as . Due to the superheating and supercooling phenomena (it will be more obvious with increasing surface tension), the melting and freezing temperatures may be not just at very small neighbors of Tc, and might be not symmetric about Tc, thus an uncertainty occurs if simply using these direct simulations in canonical ensembles to determine the coexistence temperature.

In GCE, the whole S′(E) curve can be constructed, thus we can integrate to get S(E) and to estimate Tc from Eq. (7) in our previous works.[23–27] In this work, a much simpler approach can be applied to determine the coexistence temperature without requiring to construct the whole S′(E) curve, but a single individual simulation. In the middle of the two-phase coexistence energy region, we have

Here x is the fraction of the s phase in the coexistence state. ΔSi and ΔEi are the additional contributions of entropy and energy due to the interface between the two phases. The area of interface, A = κ(y) L2, is dependent on the fraction x and the shape of the interface. Here κ is the shape factor, y is the smaller one of x and 1 – x, and L is the length of the simulation box. For a flat plane interface, κ = 1 and independent of x, is the minimal one in all shapes of interfaces, when two phases have comparable fractions, i.e., x ∼ 0.5, then y ∼ 0.5. For example, for a spherical interface, κ(y) = (36 π)1/3 y2/3 ∼ 3, and for a cylinderical interface, κ(y) = (4 πy)1/2 ∼ 2.5. Therefore, we know that the interface is almost a flat plane in the middle of the phase-coexistence region. Thus we have that S′(E) is approximately a constant from Eq. (9), equals to 1/Tc, see Fig. 1. Therefore, we can directly get the S′(E) then the Tc via a single individual GCE simulation in the phase-coexistence states with the comparable fractions thus a flat-plane interface. This is the main idea of the presented method in this work.

Fig. 1. The gNPT simulation gives the almost constant value of the temperature in the middle of the phase-coexistence energy region as the phase transition temperature Tc = 0.619±0.003ε/kB for LJ model. Here the temperature as y axis is the calculated interior temperature of the system in the simulations, 1/S′(E), rather than an input parameter.
3. Models and methods
3.1. Models

We assess the method in determining the coexistence temperature of the Lennard–Jones fcc crystal, and that of ice Ih in the all-atom TIP4P-2005 and TIP4P-ICE water models, as well as the coarse-grained mW water.[7]

3.2. Simulation details

The simulation system is built with 32480 or 4584 water molecules for mW (approximately 400 Å × 56 Å × 43 Å) and all-atom water models (approximately 200 Å × 26 Å × 29 Å), respectively. The ice phase is generated from 20 ns equilibrated bulk ice Ih at 260 K, with x–y plane set as ice Ih basal plane, leaving the secondary prismatic plane () in contact with a bulk liquid water phase. The two-phase systems are then equilibrated for another 100 ps. The size in the x direction is about 7 times of that in the y and z directions to ensure stable interfaces and study the finite-size effect. For the LJ system, a typical coexisting system composed of 8200 fcc solid and 8000 liquid particles is built in a similar way with the 001 surface of solid as the contact plane.

All MD simulations are carried out in the GCE embedded isothermal-isobaric (NPT) ensembles (gNPT), employing LAMMPS molecular dynamics code with our added GCE module.[23] A Nose–Hoover barostat[29] with 1 ps relaxation time is used to keep the pressure at 1 atm. For all-atom water models, a time step of 2 fs is used in all simulations. The pair interaction between the oxygen atoms is truncated smoothly at 13.0 Å. Nontruncated electrostatic interactions are treated by the particle–particle particle mesh solver (pppm) with a real space cutoff of 13.0 Å and precision tolerance of 10−5. The simulation parameters in the mW water model are set normally, similar to that in reference.[7] In the LJ systems, a cutoff of 2.5σ is used to truncate the interaction.

4. Results and discussion

We first illustrate the gNPT simulations of the coexistent solid/liquid LJ system to determine the phase transition temperature and its dependence on the parameters of gNPT.

The study starts with a set of gNPT simulations of LJ particles with varied parameter H0/N to visit the states with different regions of enthalpy, corresponding to the different temperatures in pure phases and the coexistent phases with different fraction of two phases. The obtained mean enthalpy of each particle in every gNPT simulation, H/N, is from –7 to –5 (unit is the ε of the LJ potential). Here N is the number of particles. The result can be easily identified as a pure phase or a phase-coexistent state, as shown in Fig. 1. This result outlines the key idea we proposed to extract the phase coexistent temperature in an easy, fast, and accurate manner.

The coexistent temperature Tc is approximately extracted as the temperature in the middle of the coexistent region, or as the statistical average of the corresponding temperature from the platform region (–6.4 ⩽ H/N ⩽ –5.8).

We also carry out the gNPT simulations with water models. Water/ice coexisting systems are first simulated at 250 K for 100 ps to relax the interfaces and to estimate the parameter of gNPT, H0/N. For the mW model, the mean enthalpy in the coexistent region is around –10.7 (the unit is always kcal/mol in this work for water). Then we run gNPT simulations each for 20 ns with the mean enthalpy ranged from –9.9 to –11.5 are shown in Fig. 2. The result shows that the gNPT simulations locate in the phase coexistence region corresponding the enthalpy –11.3 ⩽ H/N ⩽ –10.3. Meanwhile, when H/N = –9.9,–10.1,–11.5, the gNPT simulations correspond to the NPT simulations of the bulk liquid water at 312 K, 286 K and the bulk ice at 268 K, respectively.

Fig. 2. The stable conformations in the gNPT simulations with distinct enthalpies, from all liquid water, the ice/water coexistence with different fractions of ice and water, to the complete frozen ice Ih.

It is interesting to see the temperature evolution in the gNPT simulations, as shown in Fig. 3. It takes a comparably short time (less than 2 ns) for mW water to attain equilibrium. Data between 5 ns and 60 ns is sampled for statistical average, and here we only present the former 5 ns data to show that the temperature in the 6 distinct phase coexistence systems (H ranges form –11.3 to –10.3) spontaneously evolves into the same value (275.1 K), i.e., the coexistence temperature. In Fig. 4, we compare the obtained temperature in the gNPTs with the melting point of the model 274.66± 1 K in literature[7] (color in cyan), the result shows good agreement. We can get the coexistence temperature by using any one of the gNPT in the middle region of coexistence with deviation in order of 1 K. By averaging the results in the six samples covering the water/ice coexistent system, we have Tc = 275.1 ± 0.3 K.

Fig. 3. The evolution of obtained temperature in gNPTs with distinct parameter. The initial conformation is a coexistent state of ice Ih/water. Only the first 5 ns simulation trajectories are plotted.
Fig. 4. The temperatures obtained in the gNPT simulations. Cyan region refers to 274.6 ± 1 K, the melting point of the mW water model proposed in previous study.[7] Tc = 275.1 ± 0.3 K is the average of six simulations located in the phase-coexistence region.

It has to be point out that the melting temperatures we extracted from the gNPT ensembles at different assigned enthalpy slightly differ. With four models applied in this study shown in Figs. 1 and 4–6, we obtain a subtle but affirmatively varying of the phase transition temperature in the coexistent system with different phase-to-phase fractions. This dependence should not be considered as merely a finite-size effect, since we find no strong correlation between the deviation and the applied size of the simulation systems. We look forward a more comprehensive work to elucidate this effect in the future. The deviation is small (usually smaller than 1 K), so that we simply extract one of these value as the coexistent temperature with a 1 K-level error, which already makes it become one of the most precise methods in determining the phase transition temperature nowadays.

Fig. 5. The gNPT method determined the phase transition temperature of the TIP4P-2005 water model. Cyan region refers to 252 ± 5 K from the free energy calculations, yellow region refers to 249±3 K of the direct coexistence route.[21]
Fig. 6. The gNPT method determines the phase transition temperature for the TIP4P-ICE water model. The melting point (272.2 K) proposed in previous work[30] is plotted as the dash line.

We also give the results in two all-atom water models, as shown in Figs. 5 and 6. A smaller system composed of 4584 water molecules is applied, and gNPT simulations with 40 ns at most are run.

For the TIP4P-2005 water model, the obtained coexistence temperature is 250.0 ± 0.6 K, matching with Tc(DC) = 249 ± 3 K, Tc(FE) = 252±5 K, and Tc(DC) = 249.5±0.1 K from previous researches.[5,21] The temperatures in the three simulations with H = –4.25,–4.16,–4.06 come to almost the same value within 20 ns. The further 20 ns data are applied to estimate the coexistent temperature. The TIP4P-ICE water model is known to have a better performance in generating ice/water phase diagram. Its result Tc = 270.2±0.3 K, as shown in Fig. 6, has an excellent consistency with references.[5,30] All results and comparisons with literatures are shown in Table 1.

To further explore the performance of this approach, we carry out simulations with different sizes of systems. Three gNPT simulations with 44800, 5971, 2030 water molecules are built, corresponding approximately to the sizes of 151 Å × 156 Å × 58 Å, 54 Å × 58 Å × 58 Å, 52 Å × 56 Å × 22 Å, respectively. These results are colored in red, as shown in Fig. 7. They are respectively 275.1 ± 0.3 K, 275.0 ± 0.4 K, and 275.2 ± 0.6 K, while the black one 275.1± 0.3 K is the data presented in Fig. 4. The error bar only slightly increases as the size of the simulation systems obviously decreases, and in the smallest system with 2030 molecules, the precision of the obtained melting temperature is still very good.

Fig. 7. Finite-size effect in the gNPT method determines the phase transition temperature. Data from Fig. 4 is shown in black, the other data is extracted from the simulations with 44800, 5971, and 2030 water molecules, respectively.
Table 1.

The coexistence temperature of the mW, TIP4P/2005, TIP4P/ice water and LJ models as obtained from different methods (free energy calculations, Hamiltonian Gibbs–Duhem integration, direct coexistence technique, and the current GCE (gNPT) approach).

.
5. Conclusion

It has been found that the GCE (or gNPT) approach is very applicable of sampling metastable and unstable states of thermodynamics, as we did in this work sampling phase coexistence in the ice/water and the LJ solid/liquid systems. The GCE method can be applied to obtain directly the coexistent temperature of two phases even via a single individual simulation. Good precision of the approach can be reached even in a smaller simulation system, and it is better than (at least in comparable with) the other methods in the literatures with much less calculation time. Therefore, the GCE method is a highly efficient, accurate, and easy-to-use approach in determining the phase coexistent temperature.

Reference
[1] Matsumoto M Saito S Ohmine I 2002 Nature 416 409
[2] Bai G Gao D Liu Z Zhou X Wang J 2019 Nature 576 437
[3] Russo J Romano F Tanaka H 2014 Nature Mater. 13 733
[4] Li T Donadio D Russo G Galli G 2011 Phys. Chem. Chem. Phys. 13 19807
[5] Conde M M Rovere M Gallo P 2017 J. Chem. Phys. 147 244506
[6] Gao G T Zeng X C Tanaka H 2000 J. Chem. Phys. 112 8534
[7] Molinero V Moore E B 2009 J. Phys. Chem. 113 4008
[8] Moore E B Molinero V 2011 Nature 479 506
[9] Sanz E Vega C Abascal J L F MacDowell L G 2004 J. Chem. Phys. 121 1165
[10] Sanz E Vega C Abascal J L F MacDowell L G 2004 Phys. Rev. Lett. 92 255701
[11] Smit B 1992 J. Chem. Phys. 96 8639
[12] Jorgensen W L Chandrasekhar J Madura J D Impey R W Klein M L 1983 J. Chem. Phys. 79 926
[13] Vega C Abascal J L 2011 Phys. Chem. Chem. Phys. 13 19663
[14] Kroes G J 1992 Surf. Sci. 275 365
[15] Vega C Sanz E Abascal J L F 2005 J. Chem. Phys. 122 114507
[16] Ghoufi A Goujon F Lachet V Malfreyt P 2008 J. Chem. Phys. 128 154716
[17] Vega C de Miguel E 2007 J. Chem. Phys. 126 154707
[18] Alejandre J Tildesley D J Chapela G A 1995 J. Chem. Phys. 102 4574
[19] Ladd A J C Woodcock L V 1977 Chem. Phys. Lett. 51 155
[20] Bryk T Haymet A D J 2002 J. Chem. Phys. 117 10258
[21] Conde M M Gonzalez M A Abascal J L F Vega C 2013 J. Chem. Phys. 139 154505
[22] Abascal J L F Vega C 2005 J. Chem. Phys. 123 234505
[23] Xu S Zhou X Ou-Yang Z C 2012 Commun. Comput. Phys. 12 1293
[24] Jeong S Jho Y Zhou X 2015 Sci. Rep. 5 15955
[25] Zhao L Xu S Tu Y S Zhou X 2017 Chin. Phys. 26 060202
[26] Yin L Xu S Jeong S Jho Y Wang J Zhou X 2017 Acta Phys. Sin. 66 136102 in Chinese
[27] Xu S Zhou X Jiang Y Wang Y T 2015 Sci. China Phys. Mech. 58 590501
[28] Zhang C B Ye F F Li M Zhou X 2019 Sci. China-Phys. Mech. Astron. 62 67012
[29] Hoover W 1985 Phys. Rev. 31 1695
[30] Abascal J L F Sanz E García Fernández R Vega C 2005 J. Chem. Phys. 122 234511
[31] Broughton J Q Gilmer G H 1986 J. Chem. Phys. 84 5741
[32] García Fernández R Abascal J L Vega C 2006 J. Chem. Phys. 124 144506